Quantum multimode model of elastic scattering from Bose Einstein condensates 
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Mean field approximation treats only coherent aspects of the evolution of a Bose Einstein con- 
densate. However, in many experiments some atoms scatter out of the condensate. We study an 
analytic model of two counter-propagating atomic Gaussian wavepackets incorporating dynamics of 
incoherent scattering processes. Within the model we can treat processes of elastic collision of atoms 
into the initially empty modes, and observe how, with growing occupation, the bosonic enhancement 
is slowly kicking in. A condition for bosonic enhancement effect is found in terms of relevant param- 
eters. Scattered atoms form a squeezed state that can be viewed as a multi-component condensate. 
Not only are we able to calculate the dynamics of mode occupation, but also the full statistics of 
scattered atoms. 
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A remarkably universal tool describing vast majority 
of experiments with the Bose Einstein condensates is 
the celebrated Gross-Pitaevskii equation [GPE]. It de- 
scribes a coherent evolution of the atomic mean field. 
In the Hartree interpretation, its time-dependent version 
assumes that each atom of the system undergoes identi- 
cal evolution. This is a good assumption since in typi- 
cal experiments the wave-packet of the system contains 
many thousands of particles in the same state. To use a 
term borrowed from quantum optics, the time-dependent 
GPE describes stimulated processes. In some experi- 
ments [H, however, there is a clear evidence of sponta- 
neous processes. For example, in a collision between two 
condensates, some atoms from colliding quantum matter 
droplets would inevitably scatter away from them. This 
is a loss process, which is not accounted for by the con- 
ventional GPE. Description of such phenomena calls for 
use of quantum fields instead of c-number wave- functions. 
This is not easy since, in general, field equations are non- 
linear. Instead of quantum fields, several groups used 
classical stochastic fields to imitate quantum initiation 
of spontaneous processes 0- At this point it is hard 
to access the accuracy of these methods. Solid results 
so far has only been obtained within perturbation the- 
ory ^ 1^ ' It is the purpose of this Letter to present 
the first exact nonperturbative calculation of colhsional 
losses, valid in the regime of Bose enhancement. Our 
model assumes spherical nonspreading Gaussians for the 
colliding wave-packets. No doubt it will serve as a bench- 
mark test of validity of various approximate schemes in- 
cluding classical stochastic fields. 

A system of Bosons interacting via contact potential is 
described by the Hamiltonian 



H = - d^r-9Ur,t)——^ir,t) 
J 2m 



+ ly d'r¥{r,t)¥{r,t)^{r,t)4>{r,t), (1) 

where 'l'(r) is a field operator satisfying equal time 
bosonic commutation relations, m is the atomic mass 



and g determines the strength of the inter-atomic inter- 
actions. Since the Hamiltonian |^ is of the fourth order 
in ^P, the Heiscnberg equation governing the evolution 
of the field will be nonlinear and thus, in general, ana- 
lytically and numerically untractable. However, for some 
physical systems, a Bogoliubov approximation can be ap- 
plied leading to linear Heisenberg equations. The idea 
underlying this approximation states that for some cases 
the field operator might be split into two parts and 
d. First contribution describes macroscopically occupied 
field and since its fluctuations are usually small, its oper- 
ator character might be dropped {ip becomes a c-number 
wave- function satisfying GPE). The second part S, repre- 
senting fluctuations, will require full quantum mechanical 
treatment, but as long as we neglect its back-reaction on 
"tp the evolution of S will be linear. 

In this Letter we consider a process of collision of 
two strongly occupied Bose Einstein condensates. Initial 
state of the system consists of two counter-propagating 
atomic wave-packets and the "sea" of unoccupied modes. 
For such a system the Bogoliubov approximation can be 
applied. The splitting of the bosonic field is performed 
in the following manner: 



*(r,t) = V^Q(r,t)+V-Q(r,t)+(5(r,i), 



(2) 



where the subscript ±Q denotes the mean momentum of 
the colliding condensates. Upon inserting Eq. (0) into 
the Hamiltonian ^ one obtains a collection of different 
terms. We keep only those, that lead to creation or an- 
nihilation of a pair of particles 



H = - J dV(5t(r,t) 



2m 



(3) 



One can argue that such an approximation gives correct 
results if and only if the kinetic energy associated with 
the center-of-mass motion is much larger than the inter- 
action energy per particle, ti^Q^ /{2m) 3> gn, where n is 



2 



the average density of the particles in the condensates. 
Numerical proof of the above statement in the simplest 
case of two plane matter waves was given in Q • This con- 
dition is readily fulfilled in current experiments Q, |Bi 
and all the results below are obtained in this regime. 

In order to further simplify the dynamics we compare 
three characteristic timescales that appear in the prob- 
lem; the collisional time, tc = (ma) /{HQ), the time it 
takes for each wave-packet to pass through its collid- 
ing partner, the linear dispersion time, t^u — ma^ /% 
, characteristic time of the spread of the wave-packet 
due to kinetic energy term (neglecting the nonlinearity) , 
and nonlinear dispersion time, t^vD = \/ ir'l'^mcr' j gN ^ 
time of ballistic expansion in Thomas Fermi approxima- 
tion 1^. Here each of the wave-packets has the radius 
of a and contains N /2 atoms. The dynamics of our sys- 
tem depends on the relations between timescales defined 
above. Hence we introduce dimensionless parameters: 
tLo/tc = /3 and {Ild /tNoY' = a. When the num- 
ber of elastically scattered atoms is small in comparison 
with the total number of atoms in both wave-packets and 
both linear and nonlinear dispersion timescales are much 
longer than the collisional time {{tLo/tc) = /? ^ 1 and 
{tMo/tcY = ^ 1), we can neglect the change of 

shape of the macroscopically occupied functions ?/'Q(r, t) 
during the collision. In our model we use spherically 
symmetric Gaussian wave-functions 
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where r — [xi,X2,x^). In the dimensionless units, 
[t = t/tc and Xi = Xi/a, for i = 1,2,3), the Heisenberg 
evolution equation of the field operator 5 = jexp (i/3i/2) 
can be obtained upon substituting ^ into Q 



iPdt6{r,t) ^ -- {A + P^) Sir,t)+ae-'-'-* SHr,t). (5) 



The above equation has spherical symmetry! Hence, we 
decompose S into the basis of spherical harmonics 

<5(r,t) = ^ i?„,i(r)yi™(0,0Xi,„.(t), (6) 



where an,i,m are annihilation operators for a particle in 
the mode described by 71,1,111, quantum numbers. There 
is still a freedom of choice with regards to the set of 
orthogonal functions i?„,;(r). As we shall see below a 
good candidate is a set of eigenfunctions of spherically 
symmetric harmonic oscillator. 



i?n,i(r) 



Ln 



(7) 



where L'jt ^ (x) is the associated Laguerre polynomial 
and ao, a harmonic oscillator length, is an auxiliary free 
parameter that can be chosen to minimize the computa- 
tional effort. The evolution of an,i,m{t) is described by 

- _ En,i - . in" -L 



2(3 
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where coefficients Dn.i = v?'^("^ + ^ + V2)/(2/3ao)' 
En J = (2n + l + 3/2)/ag and 



Cn,n',i = / Hdri?„.;(r)exp(-r'^)i?„,;(r) 
Jo 



T{n + l + ^)r{n' + 1 



n-\-n' 



F 



(^-n,-n',l+ 



(9) 



Here F{a,b,c,x) is a hypergeometric function 0. No- 
tice that all coupling coefficients are calculated analyti- 
cally and the dn,i,m operators for different I and m are 
decoupled. Moreover, equations © do not depend on 
quantum number m. With all these simplifications the 
linear system of equations ||SJ| can be solved numerically. 

The solution of the set of dynamical equations |H1 for 
<in,i,m contains the full information about the considered 
quantum system. In particular, we can reconstruct the 
operator 5{v,t), using decomposition defined in Eq. ((HJ. 
The most straightforward observable quantity, the num- 
ber of elastically scattered atoms as a function of time 
can be expressed in terms of the trace of the density ma- 
trix 



S{t) = I d^r(SHr,t)S{r,t)) = 

Y.Y.^2l + l){al,^„^{t)d,aM, (10) 



n=0 1=0 



where {21 + 1) accounts for the degeneracy of Eq. © 
with regards to the quantum number m . In the limit 
where a/ (3 is small (notice that in Eq. this coefficient 
multiplies the source term) , S it) can be evaluated in the 
first order perturbation approximation giving Q 



ttq; 



(11) 



The same result is obtained using imaginary scattering 
length method 0]. Quality of this approximation is il- 
lustrated in Fig^ 

The bonus of having solved the full set of operator 
equations is that calculating full density matrix of the 
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FIG. 1: Number of scattered atoms versus time in pertur- 
bative regime; dashed line - analytical result given by 1)11^ . 
Solid line - numerical result obtained from our model (using 
Eq. IIUH 'I for a = 20 and /3 = 60. The inset shows the time 
evolution of the largest eigenvalue of the density matrix. 



FIG. 2; Number of scattered atoms versus time in non- 
perturbative regime where the bosonic enhancement occurs. 
Dashed line - analytical result given by CIJ- Solid line - nu- 
merical result obtained from IIOII . Parameters are: a = 160, 
/3 = 40. The inset shows the time evolution of the biggest 
eigenvalue of the density matrix. 



system of scattered atoms (p(r,r',t) = {S^r,t)5{r',t))) 
or even higher order correlation functions is just as easy 
as finding S{t). In the basis ((HJ, due to the decoupling 
property, density matrix can be written as a direct prod- 
uct oi pn,n',i,m = (oj, ; „ (i)a„' (t)) matrices, for differ- 
ent I and m. In the inset of FigQlwe present the time 
evolution of the largest of the eigenvalues of the density 
matrix /9(r, r',i). Due to the normalization of the den- 
sity matrix, J2i Kit) — S{t), where A,; are the eigenvalues 
of the density matrix, the inset of Fig^ shows that for 
a = 20, /3 = 60 there is much less than one particle even 
in the mostly populated eigenmode. 

Figure 121 shows analogous comparison between pertur- 
bative solution Ullfl and formula Hl()|l in the regime of 
parameters where the perturbation theory is expected to 
fail (the criterion for bosonic enhancement is a//? > 1 
p^2]). The figure shows that until some critical time, ap- 
proximately equal to 0.2 tc both the perturbative and 
full solutions agree very well. At this critical time the 
formula l|l()(l exceeds the perturbative solution and the 
difference between curves rapidly grows in time. At the 
same time the biggest eigenvalue of the density matrix 
of the system reaches one, which means that there is one 
particle in the mostly populated eigenmode. This ob- 
servation gives explanation to the growing discrepancy 
between two curves shown in Fig|21 Once approximately 
one atom is scattered into one of the eigenmodes of the 
density matrix the probability of scattering another atom 
into this mode grows rapidly. This is due to bosonic 
statistics of the scattered atoms and is called bosonic en- 
hancement effect. 

An interesting information about the system might 
be obtained upon analyzing the largest eigenvalues of 



Pn.n'.i,m for cach quantum number I. Figure|21juxtaposes 
these eigenvalues as a function of I, for the case with 
bosonic enhancement. The plot shows that the density 
matrix has several eigenvalues of the same order. Such a 
system is similar to quasi-condensate, in contrast to the 
commonly used definition of the condensate as described 
by a density matrix having one dominant eigenvalue |lll | . 
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FIG. 3: The biggest eigenvalue of the density matrix for dif- 
ferent / for a — 160 and (3 = 40, at time t — 2tc- Several 
eigenvalues of the same order indicate the presence of the 
quasi-condensate. 

From the experimental point of view, coherence prop- 
erties of the scattered atoms are of great importance. 
These properties are best characterized by the correla- 
tion functions. In particular, the first and second order 
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FIG. 4: First and second order correlation functions in mo- 
mentum space t;i(k, k') and (;2(k,k') for |k| — |k'| = Q, as a 
function of relative azimuthal angle 6 at t = 2tc for a — 160 
and /3 = 40. 



correlation functions can be measured in experiment. In 
one of the most commonly used method, time-of-flight 
measurement, the momentum distribution of the system 
is obtained. Thus here we calculate the first and second 
order correlation functions in momentum space using 



5i(k,k',t) 



for the former, and 



(,5t(k, t)S{k',t)) 



=(12) 



52(k,k',i) = 



{S^k,t)S^k',t)S{k',t)S{k, t)) 

{s^k,t)s{k,t)){SHk\t)S{k',t)) 



(13) 



for the latter. Due to spherical symmetry of Heisenberg 
equation for 6, the momentum density {6'^ {k,t)S{k,t)) 
is spherically symmetric as well. Moreover, since the 
Hamiltonian (|3J) is quadratic in S and the initial state 
is a vacuum state, than, in Schrodinger picture, at any 
later time t the state of scattered atoms is a multimode 
squeezed state According to general properties of 

multimode squeezed states, the n-th order correlation 
function g„(k, k') for k = k' is equal to n\. It is con- 
firmed by our numerical results. The solid line in Fig0| 
shows the first order correlation function H12(l plotted for 
fixed length of the k and k' vectors (|k| = |k'| = Q) as a 
function of relative angle 9. As expected, for 6* = the 
condition, gi{k,k) = 1 is satisfied. Also, the limited co- 
herence angle, due to spontaneous initiation of scattering 
process is clearly visible. The dashed line Fig0]shows the 
second order correlation function (|13l) . Once again, a pre- 
diction (72(k, k) = 2 is met. As Fig0]shows, the 52 func- 
tion reveals strong correlation between atoms scattered 
in direction k and — k which corresponds to relative an- 
gle 9 — 180°. This is an intuitive result, since atoms get 



scattered in pairs in such a way that the momentum and 
energy conservation laws are satisfied. Finally, the width 
of the correlation peak of §2 in the forward direction in 
the perturbative regime scales as which is propor- 
tional to the size of colliding wave-packets This 
is in analogy to Hanburry-Brown and Twiss method of 
estimating sizes of distant stars by measuring intensity- 
intensity correlation function |l.3l | and relating density- 
density correlation of 7r-mesons to the size of fireball in 
high energy collision of hadrons 0| . 

In conclusion, upon analyzing the quantum model of 
two counter-propagating atomic Gaussian wave-packets 
we get a deeper insight into processes of elastic collision 
losses of atoms and are able to study the transition from 
spontaneous regime to the bosonic enhancement. Scat- 
tered atoms form a squeezed state that can be viewed 
as a multi-component condensate. Within this model in 
principle all order correlation functions are accessible and 
hence it has a high predictive power. 
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